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ABSTRACT 

/ 

The  design  of  modern  manufacturing  systems  presents  a  number  of 
challenges.  In  particular,  the  stochastic  nature  of  machine  failures  in 
combination  with  the  large  number  of  decision  variables  makes  optimization  of 
such  systems  difficult,  'j  In^  this  paper ^we  present^ two  new  approaches  to 
optimization  of  the  complex  stochastic  systems  that  arise  in  a  manufacturing 
context;  both  are  Monte  Carlo  simulation-oriented,  and  are  therefore  broadly 
applicable.  The  first  technique  involves  using  a  likelihood  ratio  gradient 
estimate  to  drive  a  Robbins-Monro  algorithm,  and  is  relevant  to  problems  in 
which  the  decision  variables  are  continuous.  The  second  idea  employs  homotopy 
methods  to  follow  an  "optimal  path"  in  decision  variable  space,  and  can  be 
used  for  both  discrete  and  continuous  optimization. 
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MONTE  CARLO  OPTIMIZATION  OF  STOCHASTIC  SYSTEMS:  TWO  NEW  APPROACHES 


Peter  W. 


Glynn1  and  Jerry  L. 


Sanders 


1 .  INTRODUCTION 

The  problem  of  deriving  an  optimal  system  design  for  a  modern  automatic 
assembly  or  manufacturing  system  is  a  formidable  task.  A  machine  may  involve  up 
to  100  separate  assembly  or  fabrication  stations.  For  each  station,  there  are 
design  variables  which  affect  the  corresponding  cycle  time  and  buffer  area.  In 
addition,  there  are  questions  of  station  sequencing,  the  number  of  pallets  to  be 
placed  on  the  line,  whether  or  not  parallel  processing  is  to  be  used  and  if  so 
of  what  type.  Furthermore,  when  inspection  and  repair  are  required,  one  must 
decide  whether  or  not  to  use  a  multiple  loop  system  or  to  include  the  functions 
in  a  single  loop  design.  Finally,  there  are  a  variety  of  gating  and  control 
rules  that  must  be  selected  for  implementation  by  the  programmable  controllers 
that  operate  various  sectors  of  the  machine.  In  total,  there  may  be  as  many  as 
200  to  300  decision  variables  that  will  affect  the  machine's  ultimate 
productivity.  In  addition,  the  complexity  of  the  problem  is  greatly  increased 
by  the  stochastic  nature  of  the  station  failures  in  such  systems. 

Optimizing  the  design  of  these  complex  machines  is  therefore  a  major  focus 
of  current  research  in  the  area.  In  this  paper,  we  will  present  two  new 
approaches  for  optimizing  the  type  of  complex  stochastic  systems  which  arise  in 
manufacturing.  Roth  methods  are  Monte  Carlo  simulation-oriented,  and  are 
therefore  applicable  to  an  enormous  variety  of  different  problems.  In  Section 
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2,  we  describe,  in  a  precise  sense,  the  stochastic  optimization  problem  to  be 
considered  here.  Sections  3  through  5  discuss  our  new  techniques,  while  Section 
6  summarizes  our  results. 

2.  FORMULATION  OF  THE  STOCHASTIC  OPTIMIZATION  PROBLFM 

An  essential  requirement  for  the  development  of  stochastic  optimization 
algorithms  is  a  sound  probabilistic  framework  for  their  analysis.  The  following 
problem  will  serve  to  motivate  our  discussion. 

(2.1)  EXAMPLE.  Consider  an  M/M/1/*  queue  (see  [3],  p.  273)  in  which 
customers  arrive  at  rate  X .  The  goal  of  the  engineer  is  to  determine  an 
optimal  service  rate  u • 

We  formalize  the  stochastic  optimization  problem  as  follows.  Let  0  € 
represent  the  decision  variable,  and  let  KC  rf5  correspond  to  feasible  choices 
of  the  decision  variable  8 .  As  is  standard  in  the  analysis  of  stochastic 
systems,  ft  shall  represent  the  set  of  possible  sample  outcomes  w  of  the 
stochastic  system  and  F  corresponds  to  a  suitable  ff -field  of  subsets  (events) 
contained  in  ft. 

(2.1)  EXAMPLE  (continued).  Here  d  *  1,  8  =  u ,  and  X  =  (0,«).  The  sample 
space  ft  consists  of  the  set  of  right-continuous  functions  w  *  uj  (•  )  with  left 
limits,  taking  values  in  Z+  =*{  0, 1 ;  the  function  io(»)  should  be  viewed 
as  the  queue-length  process.  Specifically,  ft  *  D[0,<»)  and  F  is  the 
associated  Borel  o-field  (see  [4],  p.  116-127). 

The  concept  of  optimality  requires  an  objective  function  of  some  kind. 

Given  the  outcome  u>  and  decision  parameter  0,  let  Z(8)  =  Z(0,oj)  be  the 
"cost"  incurred  by  the  system. 

(2.1)  EXAMPLE  (continued).  Suppose  that  the  system  pays  out  one  dollar  per 
unit  of  time  that  each  customer  is  in  the  system,  and  that  the  server  is  paid 


bp  dollars  per  unit  of  time.  The  long-run  cost  of  operating  the  queue  is  then 
given  by 


Z(9,w) 


11m  1 

t*« 


—  /  id  (s)ds  +  b0 
*  0 


The  distribution  on  ft  is  clearly  affected  by  the  choice  of  9 . 
Specifically,  the  probability  measure  P0  on  (£),F)  is  allowed  to  depend  on 
9. 

(2.1)  EXAMPLE  (continued).  The  probability  P0  is  that  associated  with  a 
continuous-time  Markov  chain  on  Z*  having  initial  distribution  p (9 )  and 
generator  Q(9),  where 

X  ,  j-i  +  1,  i>0 

-x  -  e,  j-i,  i  >  o 

-9  ,  j  -  i  -  1,  i  >  1 

0  ,  else  . 

The  goal  of  the  engineer  is  to  minimize,  over  the  decision  variable  9 , 
the  expected  cost  of  running  the  system.  Mathematically,  the  stochastic 
optimization  problem  is: 

rain  r(9  )  (2.2) 

9cX 


0±j(9) 


where  r(9)  -  E0Z(9)  =  /  Z(8  ,u)PQ  (du»)  . 

n 

(2.1)  EXAMPLE  (continued).  For  the  M/M/1/*  queue,  it  is  easy  to  calculate 
that  r(9)  -  X ( 9  -  X)”1  +  b9.  The  minimizer  is  then  given  by  9*  -  X  +  (X/b)1/2. 

We  shall  say  that  (2.2)  is  a  continuous  parameter  stochastic  optimization 
problem  if  X  is  an  open  set,  and  discrete  parameter  if  X  is  a  discrete  set. 


3- 


3.  CONTINUOUS  PARAMETER  MONTE  CARLO  OPTIMIZATION 


For  the  next  two  sections,  we  specialize  (for  the  sake  of  simplicity)  to 
the  case  where  the  decision  variable  is  a  scalar  (i.e.  d  «*  1).  In  the 
continuous  parameter  setting,  the  optimization  problem  (2.2)  can  be  analyzed  in 
terms  of  finding  the  root  6*  to  the  equation  r'(0)  “  0  (provided,  of  course, 
that  r  is  differentiable  with  a  minimizer  on  K).  One  iterative  method  for 
accomplishing  this  (when  r  is  convex)  then  takes  the  form 

-  9k  -  »kr,<9k>  <3-’> 

where  is  a  suitably  chosen  sequence  of  constants  tending  to  zero.  A 

difficulty  in  implementing  (3.1)  for  the  complex  stochastic  systems  under 
consideration  here  is  that  r' (0 )  is  not  derivable  in  closed  form.  In  fact, 
r(6)  is,  in  general,  not  known  and  must  itself  be  estimated  via  Monte  Carlo 
techniques.  Such  algorithms  are  well-known  and  easily  constructed;  on  the  other 
hand,  Monte  Carlo  procedures  for  estimating  the  derivative  r' (0 )  are  not 
widely  available. 

As  a  consequence,  it  seems  (at  first)  reasonable  to  limit  oneself  to 
optimization  algorithms  in  which  only  estimators  for  r(0)  are  required  (as 
opposed  to  the  derivative  r'(0)).  The  Kiefer-Wolfowitz  (KW)  algorithm  [6] 
circumvents  the  need  to  estimate  r'(0)  by  instead  estimating  a  central 


difference  approximation  to  r'(9)«  Specifically,  let 


9k  -  »k hr-  <V0k  +  ck>  -  V0k  -  ck>>] 


*l2<^  '  k'  k 


where  +  c^)  and  Y^(0k  -  c^)  are  conditionally  independent  given 

O.,,...,©^  and  E(Yk(©k  ±  ck)  !©  ...  ,0k)  -  r(©k  t  ck**ak  anfl  ck  are 
deterministic  constants  tending  to  zero).  The  convergence  rate  of  @n  to  9 
has  been  extensively  studied  (see,  for  example,  (9])  and  the  resulting  analysis 
has  shown  that  convergence  (except  In  rare  cases)  is  slower  than  0(n  '  ). 


Since  the  typical  convergence  rate  for  Monte  Carlo  algorithms  is  0(n”1//2), 
this  is  somewhat  unsatisfactory. 

Suppose,  however,  that  Monte  Carlo  estimates  for  r' (0 >  are  available. 
The  Robbins-Monro  (RM)  algorithm  [7]  is  the  direct  stochastic  analog  of  (3.1), 
namely 


e. 


’k+1  “  6k  "  akYk*®k*  (3.2) 

where  1 8 ^ , . . . ,0^)  -  r(0^).  It  turns  out  that  (3.2)  enjoys  a 

convergence  rate  of  0(n~^2)  (see  p.  381  of  [9])}  furthermore,  the  improved 


convergence  rate  of  RM  is  manifested  empirically.  This  discussion  suggests  that 
one  should  focus  attention  on  optimization  algorithms  that  iterate  on  the  basis 
of  Monte  Carlo  sampling  of  the  derivative  (as  in  RM),  as  opposed  to  the  function 
value  itself  (as  in  KW). 

To  accomplish  this  goal,  one  therefore  requires  Monte  Carlo  estimates  of 
r'(0).  Recall  that  r(0)  “  EqZ(9).  If  the  expectation  did  not  depend  on  0 
(i.e.  if  r(8)  *  EW(0)),  then  (assuming  the  interchange  of  derivative  and 
expectation  was  valid)  it  would  follow  that  r'(0)  *  EW'(0).  Thus,  by  sampling 
W*  (0  ) ,  one  would  have  a  Monte  Carlo  estimate  for  r'(0). 

To  eliminate  dependence  of  the  expectation  on  8 ,  two  basic  approaches  are 
available.  Both  are  most  easily  illustrated  when  8  »  R1,  In  which  case  (•  ) 
is  determined  by  a  distribution  function  P(0,*).  The  first  idea  Involves 
observing  that 


r (0  )  -  /  Z(0 ,x)F(0 ,dx> 


Ace1 
N?T  ’ 


i ,  :  ,T 


jJ  ■ 
IJ- 


1 


-  /  Z(0,F-1(0,x))dx 

0 


Hy 

Dl 


EZ (0  ,F_19  ,0) ) 


i  Av  i ' 

f - - 

i 


m 
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where  F"1^,*)  is  the  inverse  distribution  function  defined  by  F-1(9,x)  » 
sup{y:F(0,y)  <  x} ,  and  U  is  a  uniform  (0,1)  r.v.  Thus,  if  Z(»  ,x)  and 
F“1(*,x)  are  differentiable,  it  follows  that  r(0 )  -  EW(9)  where  W' (0  ) 
exists.  This  idea  of  representing  the  response  in  terms  of  a  single  uniform 
r.v.  is  basically  the  method  of  common  random  numbers  (see  [8]).  It  turns  out 
that  this  approach  applies  also  to  discrete-event  systems;  the  estimator  W'(0) 
developed  in  the  discrete-event  context  is  the  perturbation  analysis  estimator 
due  to  Ho  and  Suri  (see,  for  example,  [10]). 

For  the  second  approach,  assume  that  F(9 ,• )  has  a  density  f(9  ,*  ).  Let 
g(* )  be  the  density  of  any  r.v.  X  for  which  the  density  is  positive 
everywhere  (e.g.  X  ~  H(0,1>).  Note  that 

00 

r (0 )  -  /  Z (0 ,x)f (6 ,x)dx  (3.3) 

“  /  Z(9,x)  9<*>dx 

-  E(z(0,X)f (0,X)/g(X))  . 

If  Z(*,x)  and  f(*,x)  are  differentiable,  then  clearly  r(0 )  can  be 
represented  as  EW(0)  where  W' (0 )  exists.  The  idea  described  above  is  an 
application  of  importance  sampling  [8]  to  the  derivative  estimation  problem. 

For  those  familiar  with  statistics,  one  may  interpret  the  function  f(8,X)/g(X) 
as  a  likelihood  ratio;  as  a  consequence,  we  refer  to  Monte  Carlo  gradient 
estimators  based  on  (3.3)  as  likelihood  ratio  gradient  estimators. 


4.  LIKELIHOOD  RATIO  GRADIENT  ESTIMATORS 

We  now  wish  to  describe  the  likelihood  ratio  gradient  estimation  algorithm, 
as  applied  in  the  discrete-time  countable  state  Markov  chain  context.  Assume 
the  state  space  S  of  the  chain  is  a  subset  of  Z+.  Here,  we  let  the  sample 


-6- 


space  ft  ■  {«  ■  (Wq,(i>  ,,...)  :u>1  e  S}  .  Analogous  to  the  situation  in  Example  2.1, 
a  probability  measure  Pg  is  induced  on  ft  by  specifying  an  initial 
distribution  y(0)  -  (u(0,i)si  e  S)  and  a  transition  matrix 
P(0)  »  (P(0,i, ' )si, j  €  S).  The  probability  Pg  is  then  defined  by 

n-1 

Pe{x0  "  i0'-**'Xn  “  V  “  l,(e'10>  *  n  p<0'1k'ik+1)  {4,1) 

k-0 

where  X1(u)>  is  the  coordinate  r.v.  given  by  X1(o)>  -  o>i.  We  assume  throughout 
our  discussion  that  y(0)  and  P(0 )  are  ( element-by-e lament )  continuously 
differentiable  in  0  in  a  neighborhood  of  0g,  and  that  A (0 )  ” 

{ (i, j):P(i, j,0)  >  0}  is  independent  of  0  in  a  neighborhood  of  0  0« 

The  "cost"  r.v.  Z(0,u)  will  initially  be  assumed  to  be  of  the  form 
z (0  ,(i) )  -  g(0,xo(w ),..., Xt(0)),T(U)))  where  T  is  a  stopping  time  for  X  (this 
basically  means  that  the  event  {T  ■  n}  is  determined  by  the  history  of  X  up 
to  time  n,  for  all  n  >  0)  see  (2),  p.  95)  and  g(* ,i0, . . . ,in,n)  is 
continuously  differentiable. 

(4.2)  EXAMPLE.  If  Z(6,u)  -  bXn(w),  the  cost  is  proportional  to  the  state  of 

the  chain  at  time  n. 

(4.3)  EXAMPLE.  If  Z(0,<o)  -  bT(A,w)  where  T(A)  -  inf{n  >  0:Xn  e  A},  the 

cost  is  prc  tional  to  the  firBt  hitting  time  of  A.  This  is  useful  when  one 
wants  to  determine  a  0  under  which  the  chain  will  enter  A  as  quickly  as 
possible.  In  particular,  if  A  corresponds  to  states  of  a  manufacturing  system 
in  which  the  system  is  operating  efficiently,  this  may  be  of  interest  from  a 
real-time  control  standpoint. 

Our  goal  is  to  obtain  Monte  Carlo  estimators  for  r'(0o>  •  where  r(0  )  - 
Eg Z ( 0 ) •  Now,  it  is  easily  seen  that 

-  a,  z(0  >l(0  ) 


Eg  Z(0  ) 


(4.4) 


where  L(9)  =  u(6,x0)  •  H  L(0,Xk>Xk+1)  and  L(0 ,1, j )  «  P(9 ,i, j )/P(9 Q,i, j ) . 

k=0 

(Decompose  both  sides  of  (4.4)  over  the  possible  sample  outcomes  of  (X^,..., 
and  use  (4.1).)  Thus,  r(9)  *  EU  W(9)  where 

W'(90)  “  Z* (9  0)  +  Z(0 q)L* (9 0) 

T-1 

and  L'(9g)  *  I  L' (90,Xk,Xk+i)  (observe  that  L(0g,i,j)  “1).  We  now 
k=0 

illustrate  the  method  with  an  example. 

(4.5)  EXAMPLE.  Suppose  that  we  wish  to  minimize  the  expected  time  to  empty 
an  M/M/ 1 /*  queue  qiven  that  the  oueue  initially  contains  m  customers.  Th 
queue  is  charged  b9  dollars  per  unit  time  to  run  the  server  at  rate  0 .  Th 
problem  can  be  posed  in  terms  of  the  associated  embedded  discrete-time  Markov 
chain  X  *  {x^n  >  0} .  In  particular,  let  y(0,i)  -  6 im, 

A/(A  +  9),  j-i+1,  i>1 

9/(A  +0),  j  *  i  -  1,  i  >  1 

P(9,i,j)  -  • 

1  ,  i  -  0,  j  -  1 

0  ,  else  , 

Z(8 )  -  T(1  +  b9 )/(A  +  9 ) 

where  T  *  inf{n  >  0:Xn  -  0} .  To  generate  WJ(9g): 

1)  Simulate  the  chain  X  up  to  time  T,  under  initial  distribution  y  (9  Q) 
and  transition  matrix  P(0g). 

2)  Calculate 


W}(90)  -  T 


(bA  -  1)  T(  1  +  b6) 


T-1 

l  L' (0  0,Xk,Xk+1) 
k-0 


where  L' (0 Q,i ,i  +  1 )  -  -(A  +  0 0)_1  and  L' (0 Q,i,i  -  1 )  -  A/9 Q(A  +  9 Q) . 


The  sampling  technique  described  above  could  now  be  used  at  each  step  of 
the  RM  algorithm  given  by  (3.2),  in  order  to  minimize  r(0). 


In  many  application  settings,  the  goal  is  to  minimize  r(0),  where  Z(0  ) 


is  a  steady-state  limit  (see  Example  2.1).  In  general,  the  steady-state  depends 


on  the  infinite  history  of  the  process  and  therefore  Z(0 )  cannot  be  chosen  to 


depend  on  X  only  up  to  the  stopping  time  T  (in  other  words,  Z(0 )  cannot  he 


represented  as  g{0 ,Xg, . . . ,Xtp,T) ) .  However,  in  the  current  context,  we  are  in 


luck  since  the  regenerative  structure  of  Markov  chains  allows  one  to  reduce  the 


analysis  from  the  infinite  horizon  to  that  of  a  single  regenerative  cycle. 


Specifically,  suppose 


z<8>  f(9'V  • 

k*0 


Regenerative  process  theory  [3]  shows  that  if  X  is  an  irreducible  positive 


recurrent  Markov  chain. 


Z(0 )  -  u(0 )/4(0 )  a.s. 


where  u(0)  ■  Eg(  T  f(0,X^))  and  l (0 )  -  f^T.  (Pgt*)  is  the  probability  on 
k-0 

ft  under  which  u  (0  '*>  -  « oi  and  the  transition  matrix  P(0 )  is  usedi 


T  »  inf{n  >  IjX  ■  0}  is  then  a  regenerative  stopping  time  for  x).  Hence, 


r (0 )  -u(0)/*(0>  and  r'(0)  -  [u'  (0)1(9)  -Jt'(0)u(0)]  •  £  (0  )  .  Thequantities 


£(0)  and  u(0)  can  be  estimated  using  standard  Monte  Carlo  methods.  The 


estimation  of  u'(0)  and  l' (0 )  is  easily  accomplished  using  the  analysis 


discussed  above,  since  u(0)  and  t(0)  are  the  expectations  of  r.v.'s  which 


depend  on  X  only  up  to  a  stopping  time  T.  Furthermore,  the  estimator  obtained 


in  this  way  is  suitable  for  evaluation  at  each  step  of  a  RM  algorithm  for 


minimizing  the  steady-state  cost  r(0 )  over  0 . 


Further  work  (in  particular,  extensions  to  more  general  processes)  may  be 


found  in  [5] . 


m 


v  v  v  • ,  .< 
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5.  A  HOMOTOPY  METHOD  FOR  STOCHASTIC  OPTIMIZATION 

In  this  section,  we  shall  briefly  describe  a  new  Monte  Carlo  algorithm 
which  holds  considerable  promise  as  a  technique  for  optimizing  stochastic 
systems  in  both  the  discrete  and  continuous  context.  The  idea  focuses  on  the 
fact  that  in  virtually  all  stochastic  optimization  problems,  the  system  contains 
both  controllable  and  uncontrollable  parameters. 

(2.1)  EXAMPLE  (continued).  In  the  M/M/1/“  queue,  the  system  depends  on  both 
X  (uncontrollable  parameter),  and  u  (the  controllable  decision  parameter). 

Returning  to  our  formulation  of  Section  2,  this  suggests  that  the 
probability  measure  P  on  ft  depends  both  on  the  controllable  decision 
variable  0  and  a  vector  of  uncontrollable  parameters,  say  a  e  R®,  so  that 
P  *  P(9,a).  Our  goal  is  to  find  0*  to  minimize  r(0,ao)  over  9  c  K  (a  q  is 
a  given  prescribed  setting  for  the  uncontrollable  factors),  where 

r  (0  ,a  )  *  f  Z(9  ,w)P(0  ,a  ,dw)  . 

a 

Our  basic  idea  is  that  there  may  be  a  region  T  in  the  a -space  Rf”  in  which 
the  optimization  problem  is  particularly  easy  to  solve  (either  analytically, 
numerically,  or  via  Monte  Carlo).  Specifically,  assume  that  one  can  find  a 
"good"  solution  0  at  the  point  a 1  e  r.  For  X  c  [0,1],  let 

a(X)  ■  (1  -  X)a^  +  X<Xq.  (^>ne  expects  that  the  minimizer  9  (X)  for  the  problem 
r(9,a(X))  will  deform  continuously  from  9  (0)  (i.e.  the  known  point  9)  into 

the  desired  minimizer  9*(1),  as  X  is  increased  from  0  to  1;  the 
parameter  X  is  called  the  homotopy  parameter.  The  minimizer  0*  =  0*(1)  is, 
of  course,  the  minimizer  of  the  original  problem. 

(5.1)  EXAMPLE.  Consider  an  asynchronous  automatic  assembly  system,  as 
indicated  in  Figure  1.  Pallets  containing  partially  completed  assemblies 
circulate  in  a  clockwise  fashion  on  the  machine.  Completed  assemblies  are 
offloaded  at  station  nj  at  the  next  station  (#1),  a  new  assembly  base  is 


loaded  onto  the  fixture  on  the  circulating  pallets  and  the  assembly  process 


begins  anew.  We  wish  to  optimise  the  expected  number  of  completed  assemblies 
per  hour,  as  a  function  of  the  number  of  buffer  units  b^  in  front  of  station 
i  and  the  number  of  pallets  p  on  the  machine  (so  6  -  (b^ , . • . ,b  ,p) ) . 
Although  the  machine  cycle  and  transport  times  are  assumed  deterministic  here, 
stochastic  effect  due  to  jamming  makes  this  problem  hard  to  solve. 


The  idea  is  to  observe  that  this  problem  is  greatly  simplified  when  the  jam 
rate  parameter  a  is  zero.  For  example,  if  the  machine  cycle  times  are  all 
identically  equal  to  c,  the  machine  is  then  synchronous.  Thus,  the  optimal 
throughput  1/c  is  obtained  by  setting  b^  ■  0  and  p  -  n.  The  homotopy 
algorithm  works  by  increasing  the  jam  rate  slowly  from  zero  to  the  "true"  jam 
rate  otg  >  0;  the  discreteness  of  the  problem  is  advantageous  here  since  for 
relatively  small  increments  in  X ,  the  optimizer  can  move  only  to  neighboring 
values  in  the  space  of  decision  variables  6 . 


This  type  of  approach  has  met  with  considerable  success  in  the 
deterministic  mathematical  programming  setting  (see  [1],  for  example)*  The 
primary  difficulty  in  applying  the  method  involves  following  the  path  6  (X  ) 
as  X  ranges  from  0  to  1.  Appropriate  path-following  algorithms  in  the 
Monte  Carlo  optimization  context  described  here  are  currently  under 
investigation • 

6.  CONCLUSION 

Two  new  approaches  to  Monte  Carlo  optimization  of  stochastic  systems  have 
been  presented  here.  The  first  is  based  on  utilization  of  likelihood  ratio 
gradient  estimators  in  RM  procedures,  whereas  the  second  uses  homotopy  methods 
to  follow  as  "optimal  path"  in  decision  variable  space.  Both  methods  are  under 


further  investigation,  at  both  a  theoretical  and  empirical  level 
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